#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>
#include <math.h>

#define WORDBYTES 4
#define WORDBITS 32
#define NAILS 0
#define BASE 2
#define ORDER 1
#define ENDIAN 0
#define BOUNDER 26   //832 bits = 26 bytes
extern int n18b, n18l, n52;

/*get the greatest common divisor of two numbers*/
size_t GreatestCommonDivisor(unsigned num1, unsigned num2)
{
  size_t lastGCD = 1, factor=2;

  for (; factor<=num2; factor++) {
    if ( num1%factor==0 && num2%factor==0 ) {
      lastGCD = factor;
    }
  }
  return lastGCD;
}

/*decide which algorithm to use*/
int decide (size_t sa, size_t sb, void *a, void *b, size_t *gcd)
{
  if (sa == 0 || sb == 0) {
    return 0;
  } else if (sa <= BOUNDER && sb <= BOUNDER) {
    return 1;
  } else if(sa != sb) {
    if ((*gcd = GreatestCommonDivisor(sa,sb))<=BOUNDER){
      return 1;
    } else {
      return 2;
    }
  } else {
    return 3;
  }
}

void A18base(void *a_address, void *b_address, void *c_address,\
             size_t a_size, size_t b_size, size_t *c_size)
{
  int i, j;
  unsigned c;
  unsigned long long p, q;

  /*cast addresses of a, b and c to the (unsigned *) type*/
  unsigned *int_a = (unsigned *) a_address;
  unsigned *int_b = (unsigned *) b_address;
  unsigned *int_c = (unsigned *) c_address;

  /*initialize c*/
  for (i=0; i<(*c_size = a_size+b_size); i++) {
    int_c[i] = 0;
  }

  /*do multiplication using algorithm 1.8*/
  for (i=a_size-1; i>=0; i--) {
    c = 0;
    for (j=b_size-1; j>=0 ; j--) {
      p = (unsigned long long)int_a[i] * int_b[j] + int_c[i+j+1] + c;
      q = p << 32;
      q = q >> 32;
      int_c[i+j+1] = q;
      c = p >> 32;
    }
    int_c[i+j+1] = c;
  }
}

void A52(void *a_address, void *b_address, void *c_address,\
         size_t a_size, size_t b_size, size_t *c_size)
{
  /*calculate the word number of a and b*/
  *c_size = a_size + b_size;

  int i, j, k;

  /*cast addresses of a and b to the (unsigned *) type*/
  unsigned *int_a = (unsigned *) a_address;
  unsigned *int_b = (unsigned *) b_address;

  unsigned *a1, *a2, *b1, *b2;
  size_t n;

  /*set the least significant n words of a to a2,
    the most significant n words to a1*/
  n = (a_size % 2!=0) ? (a_size +1)/2 : a_size/2;

  a1 = (unsigned *)malloc(n * WORDBYTES);
  a2 = (unsigned *)malloc(n * WORDBYTES);

  for (i=n-1, j=a_size-1; i>=0; i--,j--) {
    a2[i] = int_a[j];
  }
  for (i=n-1; i>=0 && j>=0; i--,j--) {
    a1[i] = int_a[j];
  }
  if (i == 0) {
    a1[0] = 0;
  }

  /*set the least significant n words of a to b2,
    the most significant n words to b1*/
  b2 = (unsigned *)malloc(n * WORDBYTES);
  b1 = (unsigned *)malloc(n * WORDBYTES);

  for (i=n-1, j=b_size-1; i>=0; i--,j--) {
    b2[i] = int_b[j];
  }
  for (i=n-1; i>=0 && j>=0; i--,j--) {
    b1[i] = int_b[j];
  }
  if (i == 0) {
    b1[0] = 0;
  }

  unsigned *t1, *t2, *ft1, *ft2;
  unsigned long long p, qq;
  int c = 0;
  size_t t1_size, t2_size;

  /*set t1=a1+a2*/
  t1 = (unsigned *)malloc((t1_size = n+1) * WORDBYTES);
  for (i=n-1, j=n; i>=0; i--,j--) {
    p = (unsigned long long)a1[i] + a2[i] + c;
    qq = p << 32;
    qq = qq >> 32;
    t1[j] = qq;
    c = p >> 32;
  }
  ft1 = t1;
  if (c != 0){
    t1[0] = c;
    c = 0;
  } else {
    t1 ++;
    t1_size --;
  }

  /*set t2=b1+b2*/
  t2 = (unsigned *)malloc((t2_size = n+1) * WORDBYTES);
  for (i=n-1, j=n; i>=0; i--,j--) {
    p = (unsigned long long)b1[i] + b2[i] + c;
    qq = p << 32;
    qq = qq >> 32;
    t2[j] = qq;
    c = p >> 32;
  }
  ft2 = t2;
  if (c != 0){
    t2[0] = c;
    c = 0;
  } else {
    t2 ++;
    t2_size --;
  }

  /*set w3=t1*t2*/
  unsigned *w3 = (unsigned *)malloc((t1_size+t2_size) * WORDBYTES);
  size_t *w3_size = (size_t *)malloc(sizeof(size_t));

  if (t1_size <= BOUNDER && t2_size <= BOUNDER) {
    A18base(t1, t2, w3, t1_size, t2_size, w3_size);
  } else if (t1_size == t2_size) {
    A52(t1, t2, w3, t1_size, t2_size, w3_size);
  } else  {
    A18base(t1, t2, w3, t1_size, t2_size, w3_size);
  }

  /*set w2=a1*b1*/
  unsigned *w2 = (unsigned *)malloc(2*n * WORDBYTES);
  size_t *w2_size = (size_t *)malloc(sizeof(size_t));

  if (n <= BOUNDER) {
    A18base(a1, b1, w2, n, n, w2_size);
  } else {
    A52(a1, b1, w2, n, n, w2_size);
  }

  /*set w4=a2*b2*/
  unsigned *w4 = (unsigned *)malloc(2*n * WORDBYTES);
  size_t *w4_size = (size_t *)malloc(sizeof(size_t));

  if (n <= BOUNDER) {
    A18base(a2, b2, w4, n, n, w4_size);
  } else {
    A52(a2, b2, w4, n, n, w4_size);
  }

  /*set c_array as w4/2^(32n), w4 as w4 mod 2^(32n)*/
  unsigned *c_array = (unsigned *)malloc((n+1) * WORDBYTES);

  for (i=0; i<=n; i++) {
    c_array[i] = 0;
  }

  for (i=*w4_size-n-1, j=n; i>=0 && j>=0; i--, j--) {
    c_array[j] = w4[i];
  }

  /*set w3=w3-w2-w4+c_array*/
  long long q;

  for (i=*w3_size-1, j=2*n-1, k=n; k>=0; i--,j--,k--) {
    q = (long long)w3[i] - w2[j] - w4[j] + c_array[k] + c;
    if (q < 0) {
      w3[i] = (unsigned)(q%(long long)pow(BASE, WORDBITS) + pow(BASE, WORDBITS));
      c = q/(long long)pow(BASE, WORDBITS) - 1;
    } else {
      qq = q << 32;
      qq = qq >> 32;
      w3[i] = qq;
      c = q >> 32;
    }
  }
  for (; j>=0; i--,j--) {
    q = (long long)w3[i] - w2[j] - w4[j] + c;
    if (q < 0) {
      w3[i] = (unsigned)(q%(long long)pow(BASE, WORDBITS) + pow(BASE, WORDBITS));
      c = q/(long long)pow(BASE, WORDBITS) - 1;
    } else {
      qq = q << 32;
      qq = qq >> 32;
      w3[i] = qq;
      c = q >> 32;
    }
  }
  while (c!=0) {
    q = (long long)w3[i] + c;
    if (q < 0) {
      w3[i--] = (unsigned)(q%(long long)pow(BASE, WORDBITS) + pow(BASE, WORDBITS));
      c = q/(long long)pow(BASE, WORDBITS) - 1;
    } else {
      qq = q << 32;
      qq = qq >> 32;
      w3[i--] = qq;
      c = q >> 32;
    }
  }

  /*set c_array as w3/2^(32n); w3 as w3 mod 2^(32n)*/
  for (i=*w3_size-n-1, j=n; j>=0 && i>=0; i--, j--) {
    c_array[j] = w3[i];
  }

  /*set w2 as w2+c_array*/
  for (i=*w2_size-1, j=n; j>=0 && i>=0; i--, j--) {
    p = (unsigned long long)w2[i] + c_array[j] + c;
    qq = p << 32;
    qq = qq >> 32;
    w2[i] = qq;
    c = p >> 32;
  }
  for (; c!=0 && i>=0; i--) {
    p = (unsigned long long)w2[i] + c;
    qq = p << 32;
    qq = qq >> 32;
    w2[i] = qq;
    c = p >> 32;
  }

  /*set c_address as w1*2^3n+w2*2^2n+w3*2^n+w4,
    where w1 is the uppermost n bits of w2*/
  unsigned *int_c = (unsigned *)c_address;

  for (i=*c_size-1, j=*w4_size-1; j>=*w4_size-n; i--, j--)
    int_c[i] = w4[j];
  for (j=*w3_size-1; j>=*w3_size-n; i--, j--)
    int_c[i] = w3[j];
  for (j=*w2_size-1; j>=0 && i>=0; i--, j--)
    int_c[i] = w2[j];

  free(a1);
  free(a2);
  free(b1);
  free(b2);
  free(ft1);
  free(ft2);
  free(w2);
  free(w3);
  free(w4);
  free(c_array);
  free(w2_size);
  free(w3_size);
  free(w4_size);
}

void A18large(void *a_address, void *b_address, void *c_address,\
              size_t a_size, size_t b_size, size_t *c_size, size_t gcd)
{
  size_t s_a, s_b;
  int i, j, k, l;
  s_a = a_size/gcd;
  s_b = b_size/gcd;

  /*cast addresses of a, b and c to the (unsigned *) type*/
  unsigned *int_a = (unsigned *) a_address;
  unsigned *int_b = (unsigned *) b_address;
  unsigned *int_c = (unsigned *) c_address;

  unsigned *a, *b, *c;
  unsigned *p = (unsigned *)malloc(2*gcd * WORDBYTES);
  unsigned *carrier = (unsigned *)malloc(gcd * WORDBYTES);
  size_t *s_p = (size_t *)malloc(sizeof(size_t));
  unsigned w_c;
  unsigned long long w_p, q;

  /*initialize c_address*/
  for (i=0; i<(*c_size = a_size+b_size); i++) {
    int_c[i] = 0;
  }

  /*do multiplication*/
  for (i=s_a-1; i>=0; i--) {
    /*initialize carrier*/
    for (k=0; k<gcd; k++) {
      carrier[k] = 0;
    }
    a = int_a+i*gcd;

    for (j=s_b-1; j>=0 ; j--) {
      b = int_b+j*gcd;
      c=int_c+(i+j+1)*gcd;
      /*calculate p*/
      //p=a*b
      A52((void *)a, (void *)b, (void *)p, gcd, gcd, s_p);

      //p=p+c+carrier
      w_c = 0;
      for (k=2*gcd-1, l=gcd-1; l>=0; k--, l--) {
        w_p = (unsigned long long)p[k] + c[l] + carrier[l] + w_c;
        q = w_p << 32;
        q = q >> 32;
        p[k] = q;
        w_c = w_p >> 32;
      }
      while (w_c!=0 && k>=0) {
        w_p = (long long)p[k] + w_c;
        q = w_p << 32;
        q = q >> 32;
        p[k--] = q;
        w_c = w_p >> 32;
      }

      /*put the lowermost gcd words of p to correct position in c*/
      for (k=2*gcd-1, l=gcd-1; l>=0; l--, k--) {
        c[l] = p[k];
      }
      /*put the uppermost gcd words of p to carrier*/
      for (l=gcd-1; l>=0; l--, k--) {
        carrier[l] = p[k];
      }
    }
    /*put carrier to c*/
    c=int_c+(i+j+1)*gcd;
    for (l=gcd-1; l>=0; l--) {
      c[l] = carrier[l];
    }
  }

  free(p);
  free(carrier);
  free(s_p);
}

void Product32(void *a, void *b, void *c, size_t sa, size_t sb, size_t *sc)
{
  size_t gcd;
  switch (decide (sa, sb, a, b, &gcd)) {
    case 0:
      *sc = 0;
      n18b++;
      break;
    case 1:
      A18base(a, b, c, sa, sb, sc);
      n18b++;
      break;
    case 2:
      if (sa >= sb) {
        A18large(a, b, c, sa, sb, sc, gcd);
      } else {
        A18large(b, a, c, sb, sa, sc, gcd);
      }
      n18l++;
      break;
    case 3:
      A52(a, b, c, sa, sb, sc);
      n52++;
      break;
  }
}
